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We explore use of the harmonic Einstein equations to numerically find stationary black holes 
where the problem is posed on an ingoing slice that extends into the interior of the black hole. 
Requiring no boundary conditions at the horizon beyond smoothness of the metric, this method 
may be applied for horizons that are not Killing. As a non-trivial illustration we find black holes 
which, via AdS-CFT, describe a time-independent CFT plasma flowing through a static spacetime 
which asymptotes to Minkowski in the flow's past and future, with a varying spatial geometry in- 
between. These are the first explicit examples of stationary black holes which do not have Killing 
horizons. When the CFT spacetime slowly varies, the CFT stress tensor derived from gravity is well 
described by viscous hydrodynamics. For fast variation it is not, and the solutions are stationary 
analogs of dynamical quenches, with the plasma being suddenly driven out of equilibrium. We find 
evidence these flows become unstable for sufficiently strong quenches, and speculate the instability 
may be turbulent. 



INTRODUCTION 

Due to the remarkable Anti de Sitter-Conformal Field 
Theory (AdS-CFT) correspondence [IHS], the behaviour 
of black holes in asymptotically AdS spacetimes is equiv- 
alent to the behaviour of hot plasma in certain strongly 
coupled CFTs. Since these black holes may have pla- 
nar horizons, they admit perturbations of arbitrarily 
long wavelength which give rise to the hydrodynamic be- 
haviour expected of the CFT plasma [4]-[8]. Perturba- 
tions on short scales correspond to microscopic plasma 
behaviour beyond hydrodynamics. As this currently can- 
not be computed directly in strongly coupled CFTs, grav- 
ity provides an entirely new computational tool, P| [TO]. 
This has been exploited for dynamical quenches where 
the CFT is abruptly perturbed pn4T5] and where the 
dual spacetime is determined by numerical dynamical 
gravity JMIS]. 

Here we study an analog of a dynamical quench, where 
the CFT state is time independent. We consider station- 
ary black holes dual to a time independent relativistic 
plasma flow through a static spacetime. This asymp- 
totes to Minkowski in the flow's past and future, but 
in-between the spatial geometry varies in the flow direc- 
tion. The flow, initially in equilibrium, is forced out of 
equilibrium in response to passing through the curved 
spacetime region, before returning to equilibrium after- 
wards. For slowly varying spacetimes (with respect to 
the length scale set by the local temperature) such flows 
are well described by hydrodynamics [23] . For quick vari- 
ation they probe behaviour beyond hydrodynamics, and 
are stationary analogs to dynamical quenches. 

These black holes are of a qualitatively new variety, 
being the first explicit examples of stationary black holes 
that do not have Killing horizons, and hence do not move 
rigidly [24^. The rigidity theorem states that if a station- 



ary horizon is compact, it is also Killing [25H27]. Our 
black holes have non-compact horizons and evade this 
theorem, and since the dual plasma flows in a direction 
which is not a symmetry, these horizons are not Killing. 
Other stationary non-Killing horizons have been conjec- 
tured in the AdS-CFT context; 'flowing funnels' [28.129] 
and 'plasma shocks' [30l[31]. So far only related solutions 
with Killing horizons have been found [321 [33] . 



HARMONIC EINSTEIN EQUATIONS 

Consider a Lorentzian stationary solution to the Ein- 
stein equations where the stationary Killing vector field 
T is globally timelike. We consider the purely gravita- 
tional case Rjj^i, = Kg^y^ although generalisation to in- 
clude matter is straightforward. We adapt coordinates, 
= (t,x*), so T = d/dt and the metric g^jj is, 

ds^ = -N{x){dt + Ai{x)dx'f + bij{x)dx'dx^ . (1) 

The spacetime is Lorentzian so det ^^^^ = — A^det bij < 0, 
and T is globally timelike so N{x) > and thus bij (x) is a 
Riemannian metric. In order to obtain a well posed prob- 
lem we must eliminate the coordinate invar iance. Instead 
of solving the Einstein equation, we solve the 'harmonic' 
or 'DeTurck' Einstein equations j34]436] , 

where, = g^^ (t^ — f ^ is constructed from a 

fixed reference connection f ^ on the manifold, which 
here we take to be the connection of a reference metric 
gj^iy. The two derivative part of these equations is gov- 
erned by the operator h^^didj^ and since bij is a Rieman- 
nian metric the harmonic Einstein equation is elliptic. 
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For suitable boundary conditions the whole system 
may be solved as a standard elliptic boundary value prob- 
lem. We want solutions with = 0, which is a coordi- 
nate gauge condition analogous to generalised harmonic 
gauge in dynamical numerical GR [37 , and must ensure 
our boundary conditions are compatible with this. In cer- 
tain cases one may prove must vanish [32]. In general, 
'soliton' solutions with 7^ may exist. However for an 
elliptic problem solutions are locally unique, and hence 
one may easily distinguish whether a solution found has 
vanishing or not. The system may be solved by relax- 
ation which is related to Ricci flow. Alternatively after 
discretization it can be solved by the Newton method 
given an initial guess. 

OLD METHOD FOR KILLING HORIZONS 

For a stationary black hole T = d/dt is no longer 
globally timelike, being spacelike inside the horizon, or 
outside if an ergoregion exists. Hence hij must become 
Lorentzian, and the problem inside the horizon and er- 
goregion naively appears hyperbolic. 

A previous method [35 focussed on retaining ellipticity 
by assuming a Killing horizon that rigidly moves. We 
assume Killing vectors Ra = d/dy^ exist which commute 
with themselves and T. For constants Q.^ we take K = 
T-\-Q^Ra to generate the Killing horizon and rigid motion 
of the spacetime. The metric can then be written as, 

ds^ = GAB{x){dy^ + A^{x)dx''){dy^ + A^{x)dx^) 

-^bab{x)dx''dx^ (3) 

with y"^ = (t^y^). Now Gab is Lorentzian outside the 
horizon (even in an ergoregion), and degenerates at the 
horizon or axes of symmetry of the Ra- Correspondingly 
bab{x) can be chosen to be Riemannian on, and in the 
exterior of, the horizon. These coordinates yield a slice 
of the spacetime that intersects the bifurcation surface 
of the Killing horizon. Since the principle part of R^ is 
governed by b^^dadb the p.d.e. system is elliptic posed on 
the base geometry with coordinates x^ - the 'orbit space'. 
The boundaries of this base are where Gab degenerates 
and smoothness determines boundary conditions there 
[35l [38| [39]. The surface gravity and (angular) veloci- 
ties of the horizon are prescribed in these boundary 
conditions. For regularity the reference metric must also 
have a Killing horizon at the same location with the same 
K: and Q^. Hence we may think of the reference metric 
as specifying these moduli of the solution. 

NEW METHOD FOR NON-KILLING HORIZONS 

Suppose we are interested in stationary black holes 
that do not have a Killing horizon. For a non-Killing 



horizon we cannot assume existence of a bifurcation sur- 
face and regular past horizon. In the new approach we 
now describe we no longer require the problem to be el- 
liptic. We take the general stationary ansatz ([T]) and 
pose the harmonic Einstein equations on an ingoing slice 
(analogous to that of Eddington-Finklestein) that inter- 
sects the future horizon and extends into the black hole 
interior. For the metric ([T]) in ingoing coordinates g^i, is 
regular at the future horizon, so det ^^j^ = —Ndet hij < 
on the future horizon and its exterior. Thus bij is elliptic 
in the exterior of the horizon or ergoregion (ie. where 
> 0) and is hyperbolic inside these (where A" < 0). 
Such a problem is analogous to mixed hyperbolic ellip- 
tic p.d.e.s in fluid dynamics. Whilst the problem will 
have hyperbolic character inside the horizon and ergore- 
gion we may still solve it using the Newton method as 
before. Interestingly the Ricci flow method appears also 
to work, but we will not explore that here. All compo- 
nents of = give non-trivial gauge conditions; the 
are associated to coordinate freedom in and to the 
freedom t ^ t -\- f {x'^) . 

An important difference to the old method is that since 
the problem is hyperbolic in the interior of the horizon, 
at the innermost points of our domain we impose only 
the harmonic Einstein equations and no boundary con- 
dition. The requirement that the metric is smooth in our 
domain is sufficient to ensure regularity of the horizon 
in ingoing coordinates. Starting from a smooth initial 
guess near a solution, then the Ricci flow and Newton 
method will preserve smoothness, since both update the 
metric using the harmonic Ricci tensor which will also be 
smooth. We implicitly assume the physically reasonable 
statement that asymptotic boundary conditions together 
with future horizon regularity define a locally unique sta- 
tionary black hole solution, up to moduli of the solution 
(such as mass). This is true for Killing horizons as can 
be seen from the elliptic nature of the p.d.e.s discussed 
earlier. Indeed the black hole uniqueness theorems show 
in many cases global uniqueness. For stationary non- 
Killing horizons we assume local uniqueness here, but 
emphasise we know of no proof. We note this is the basis 
of the fluid/gravity correspondence [3 [8]. It is the hori- 
zon rather than the ergosurface where smoothness must 
be imposed, even though the ergosurface determines the 
transition of character of the p.d.e.s. This is analogous to 
a stationary scalar field in the Kerr background, where 
one explicitly sees the scalar equation has regular singu- 
lar behaviour at the horizon, and hence it is smoothness 
there that constrains the solutions [40 . 

A second key difference with the old method is that the 
reference metric, while selecting the coordinate system a 
solution is presented in, no longer specifies the surface 
gravity or velocities of the horizon. These moduli must 
be fixed by appropriate boundary conditions. 

For a reader interested in implementing the method 
we provide a toy example in appendix A; finding 
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Schwarzschild assuming static spherical symmetry. The 
old and new methods are given and contrasted. 



HOLOGRAPHIC PLASMA QUENCHES 

An example application of the method above is to 
find stationary black holes that are locally asymptotically 
AdS^^ and are relevant in AdS-CFT to describe CFT sta- 
tionary plasma flows in a non-trivial geometry. These are 
Einstein metrics solving R^j^ = —pg^w We choose units 
so that the AdS length / = 1. These geometries have a 
conformal boundary whose conformal class we are free to 
specify and corresponds to the spacetime that the CFT 
is defined on. We choose a static metric, 



-dr 



dp^ + cF{p)dy^ 



(4) 



where cf{p) deforms the spatial geometry breaking the 
translation symmetry in p as, 
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(5) 



for constants and the geometry asymptotes to 

Minkowski for p ±oo. We take the CFT plasma to be 
stationary, homogeneous in and flowing from p = —oo 
to +00. We expect the plasma flow in the asymptotic 
Minkowski regions p ^ ±oo to become homogeneous, 
and correspondingly the dual to become a homogeneous 
black brane in these limits. However, since the spatial ge- 
ometry of the boundary metric depends non-trivially on 
the direction that the plasma flows in, the dual black hole 
must have velocity in a direction which is not associated 
to an isometry. Hence the plasma flow is inhomogeneous, 
and the dual black hole does not have a Killing horizon. 

We take the ingoing plasma to be subsonic with ve- 
locity vq < I/V2 and temperature Tq. Using the holo- 
graphic fluid/gravity correspondence [ll [8] provided the 
boundary metric gradients are sufficiently small, mean- 
ing P/Tq 0, then the plasma behaves as an ideal fluid. 
The first deviation from ideal behaviour is due to shear 
viscosity. Upon increasing (S/Tq towards unity, one ex- 
pects the derivative expansion of hydrodynamics to break 
down completely in the region where the boundary met- 
ric is highly curved, as microscopic physics is required to 
describe small scale plasma phenomena. As we shall see, 
this microscopic breakdown of viscous hydrodynamics is 
indeed captured by the dual gravity black hole, and these 
solutions represent stationary flowing plasma quenches. 

We write an ansatz for these metrics as. 



ds' 



{-Tdt^ + 2Vdtdz + 2Udtdp + Adz^ 



-VB{dp^Fdzf ^Sdy^^ (6) 



with the functions T^V^B^ F, A being smooth (or at 
least C^) in p and z. The locally AdS 4, boundary is at 



z = 0, and we impose the boundary conditions such that 

T = V = A = B = 1,U = F = 0,S = a{p) (7) 

there. We solve the Einstein equations and gauge con- 
dition = as a power series in z near the boundary 
at z = 0, and then transforming to Fefferman- Graham 
coordinates, identify the boundary metric as that in Q. 
We then extract the vev of the dual CFT stress tensor 
from the z^ terms in the expansions of T, U, . . . , A us- 
ing holographic renormalisation [41 . Details are given in 
appendix B. 

For regularity at the locally AdS boundary we require 
the reference metric to obey the same boundary condi- 
tions as the metric. We choose the reference metric to be 
a boosted homogeneous black brane but with S deformed 
to obey the boundary requirement so, 

S = a{p) ,T = l-cl {z/zof ,B = l^sl {z/zof (8) 

F = —Sr/B , A = —sl/B ^ V = Cr ^ U = SrCr {z / Zq)^ 

for constants zq and r, with Cr = coshr, Sr = sinhr. The 
horizon of the reference metric is at z = zq. We empha- 
sise that this is not a solution to the Einstein equations 
due to the non-trivial cr(p). 

We compactify p to a coordinate x G [—1,1] where 
dp = dx/{l — x'^)'^ . In particular this implies that p ~ 
l/(x — l)asx^l and similarly for x = —1. Since for 
a black brane perturbations should decay exponentially 
in p as p ^ ±00 we then expect in our coordinates (i.e. 
for our reference metric above) all the functions T, . . . , 6* 
will have all x derivatives vanishing as x ^ ±1. 

We work in the coordinate domain x G [—1,1] and 
z G [0, Zmax]' The solutions we find have a future horizon 
whose position in z is given as a function z = H{x) where 
< H{x) < Zmax so that horizon regularity is imposed by 



smoothness of T, U, . . . , A there. For z 



and —1 < 



X < —1 we impose the equations of motion as for the 
interior points. We emphasise that there is no boundary 
condition at z = Zmax- Finally at x = ±1 we impose 
Neumann boundary conditions for all the functions, as x 
derivatives should vanish there. 

For our solutions we expect two moduli which we can 
take as the surface gravity and velocity of the horizon 
in the asymptotic region p ^ — 00. In the dual picture, 
these are the temperature and velocity of the inflowing 
plasma which is in equilibrium. These moduli are not 
fixed by the reference metric in this ingoing method, and 
thus we must fix two pieces of information to specify a 
locally unique solution. This may be done in a myriad 
of ways, but we have found a numerically robust method 
is to fix a Dirichlet condition on U at z = z^ax and 
X = —1 and on T at z = z^ax and x = +1, instead of 
the Neumann conditions for these at those points. At 
these two points the value of U at x = —1 and T at 
X = are chosen to be those for the reference metric. 
Thus the two constants zq and r in the reference metric 



4 



become the two parameters of the solutions controhing 
the ingoing plasma temperature, Tq, and velocity, vq. 

We take the initial guess for the metric to be the ref- 
erence metric. We use finite differencing to obtain so- 
lutions, and discuss the tests of convergence in detail 
in appendix C, finding our code produces approximately 
fourth order convergence. For the resolutions used, up to 
70 X 280 in z and x, the maximum fractional local error in 
the Einstein equations outside the horizon, is better than 
~ 10~^. Hence these are very good numerical solutions 
and we have checked they are indeed Einstein metrics, as 
we require, rather than Ricci solitons. Convergence tests 
for the extraction of the boundary stress tensor (which 
depends on multiple derivatives of the metric functions) 
indicate it is accurate to better than percent level. 
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FROM HYDRODYNAMICS TO QUENCHES 

We now present data where the ingoing homogeneous 
plasma has subsonic velocity vq = 0.50, and temperature 
To = 0.24 in our units. Since the boundary theory is 
a CFT, any other temperature is related by an appro- 
priate scaling, and this value is taken for convenience. 
We choose the boundary metric to have a = 0.4, and 
we adjust f3 to move between a slowly or rapidly varying 
geometry. This value of a is sufficiently large that the 
boundary metric deformation from Minkowski cannot be 
described by perturbation theory. As we shall see, the 
deviation from homogeneous behaviour will correspond- 
ingly be large. With these data we find the dual gravity 
solution and from it extract the vevs of the CFT stress 
tensor components Ttt^Tfp^Tpp and Tyy. The conserva- 
tion equation together with tracelessness implies that all 
the information in the stress tensor is characterised by a 
single function of p. We choose to plot the (scale invari- 
ant) function defined by. 



1+^2 



(9) 



for < < 1 where v gives the local velocity of the 
plasma in the stationary frame. In figure [l] we plot this 
function for various P between 0.2 and 2. 

In this plot we show the same quantity for the 
fiuid/gravity viscous hydrodynamics with the same in- 
going fiow data - see appendix B for details. We see that 
for the smallest /3 = 0.2 the agreement of the gravity 
stress tensor with that of viscous hydrodynamics is good. 
The agreement becomes worse as P increases and higher 
derivative terms in the hydrodynamic expansion become 
important. For (3 o:^ 0(1) the hydrodynamic approxima- 
tion breaks down and we are in the quench regime. The 
bulk solutions remains perfectly smooth and allow us to 
compute the behaviour of this strongly coupled plasma 
fiow. The deviation from hydrodynamics becomes large; 



FIG. 1. Velocity v vs. p for the flows obtained from the 
dual gravity solution (solid lines) and by solving the equa- 
tions of fluid/gravity viscous hydrodynamics (dashed lines), 
for /3 = 0.2, 0.3, 0.5, 0.7 (top) and /3 = 1, 1.5, 2. {bottom). The 
plasma flows from left to right, starting in the same homo- 
geneous equilibrium state for the different /3. For small /3 we 
see good agreement with hydrodynamics. For large /3 we see 
strong deviations in the region < 1 where the boundary 
metric varies quickly. As /Sp ^ -\-oo the plasma equilibrates 
and recovers homogeneity. For /3 = 2. we see the flow is su- 
perluminal in the quench region. 



for /3 = 2 we find the local plasma velocity v becomes su- 
perluminal in a region where the metric is curved. This 
presumably indicates that for sufficient quench strength, 
/3, these fiows become unstable. We note that whilst the 
stress tensor is superluminal, all its components are well 
behaved - for example in appendix B we display (T^^). 

Interestingly we find the equilibrated outgoing plasma 
has a temperature, and hence entropy density, that is 
roughly independent of p. The same is true for the 
fiuid/gravity viscous hydrodynamics. One can see in fig- 
ure [l] that the outgoing velocities v are numerically close 
for the different /3, although they are not obliged to be 
by stress energy conservation. We emphasise that whilst 
the total entropy generated in these ffows is similar for 
different /3, the region where the spacetime is curved and 
hence this entropy is generated is very different, becom- 
ing small for large p. Hence for strong quenches the 
entropy density in the plasma is generated in a sudden 
non-adiabatic manner. 

The horizon, defined by the zero set of h{z,p) = 
z — H{p) is a null surface so that, g^^ d^hdyh = 0. This 
is an o.d.e. for H{p) which can be solved to find the 
horizon location. The null tangent to the horizon can 
be written, = ^ + Qh{p)R, where R has unit norm 
i?^ = 1, is tangent to the horizon and orthogonal to d/dt 
and d/dy. Then flnip) gives the local velocity of the 
horizon, and is plotted in figure [2j We note this is well 
behaved even for the fiow with P = 2 which has superlu- 



5 



a 



c\j 




FIG. 2. Velocity of the horizon Qh {top) and surface gravity 
Hi^ (bottom) as functions of p for the same flows as in Fig. Ill 
These functions explicitly depend on p as the horizon is non- 
Killing. Near the asymptotic regions p zboo both Qh and 
K become constant since our solutions approach homogeneous 
boosted black branes. 



minal boundary stress tensor. The boundary metric, and 
consequently the bulk metric, explicitly depend on p and 
so d/dp is not Killing. Thus the spacetime motion is not 
rigid, and hence the local velocity fin explicitly depends 
on p, rather than being constant. We also compute the 
surface gravity hc defined as V^(xiyX^) = —2f<ix^- Again 
this surface gravity is not constant, and is plotted in the 
same figure. It is also well behaved for /3 = 2. Further 
details of the solutions are given in appendices B and C. 



SUMMARY 

We have proposed a new numerical method to find 
stationary black hole solutions that do not have Killing 
horizons, and hence do not rigidly move, using the har- 
monic Einstein equations posed on an ingoing domain 
that pierces the future horizon. We have explicitly con- 
structed holographic duals to time independent plasma 
flows in a static geometry that interpolates smoothly 
between two asymptotic Minkowski regions. For gentle 
interpolations the plasma behaves as a viscous fluid, as 
predicted by the fluid/gravity correspondence. For sharp 
interpolations there is no hydrodynamic description, yet 
the dual black hole allows us to compute microscopic 
properties of this strongly coupled far-from-equilibrium 
plasma, such as the vev of the stress tensor. Such solu- 
tions are the stationary analog of dynamical quenches. 
Interestingly we find for a sufficiently strong quench 
the stress tensor vev may become superluminal in some 
region. This likely indicates that these flows become 
unstable for sufficient quench strength. It is possible 
this instability is turbulent in nature, in analogy with 



global AdS-Schwarzschild where there is thought to be 
a turbulent instability for superradiant solutions which 
correspondingly have superluminal dual plasma [42 . 

Comment: We understand that [43| numerically finds 
flowing funnel solutions which have non-Killing horizons. 
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Appendix A: Illustrative toy example 

In this appendix we illustrate the old method, which 
assumes a Killing horizon, and the new ingoing method 
described in this paper, using a simple toy example; 
numerically finding the Schwarzschild solution assuming 
spherical symmetry. The purpose is to contrast the two 
methods, and illustrate explicitly how to use them in as 
simple a context as possible. We hope this will be of use 
to a reader interested in actually implementing these 
methods in more complicated settings. 

The old Killing horizon method 

We write an ansatz for the black hole metric as, 



iBdr"^ S 



dn^ 



(10) 

where dQ^ is the line element on the unit round 2-sphere, 
A^B^S are smooth (at least C^) functions of a compact 
coordinate r G [0,1] and /(r) = 1 — r^. We require 
A = B = S = lditr = l giving asymptotic flatness. 
Regularity at the horizon r = implies that S must 
be smooth functions in r^, and then tz gives the surface 
gravity. 

We may discretize A^ B^ S on the interval [0, 1] and 
then require Neumann boundary conditions at r = 
for these functions. As an example, one might take the 
reference metric and the initial guess to be the above 
metric with A = B = 1 — f/2 and S = 1 (which is not 
Schwarzschild). On finding a solution, one obtains that 
odd derivatives at r = vanish. However this approach 
suggests there is a boundary at the horizon, which 
really there is not. A better way to think is solving the 
problem in the domain [—1,1] requiring smooth and 
even solutions so A{—r) = A{r) and similarly for B and 
S. Using finite difference or pseudo-spectral methods 
one may choose lattices with even numbers of points 
that avoid r = altogether. Then we have no boundary 
at the horizon and instead solve the problem on the 
complete t = slice representing the Einstein-Rosen 
bridge that intersects the bifurcation surface. 

The new ingoing non-Killing horizon method 



Take an ansatz with ingoing time 



ds^ = -Tdr 



-Vdtdz 



1 



rAdz' 



;dn^ (11) 



for a compact coordinate z e [0,1] with T, A, S suffi- 
ciently smooth (at least C^) functions of z. We impose 
asymptotic flatness disT = V = S = l and A = at 
z = 0. For monotonic S the horizon occurs at T = and 
provided V is non-zero and the functions are smooth in 
z there then the horizon will be regular. Now 2: = 1 is 
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not regarded as a boundary, and the equations ([2| are 
imposed there as in the interior of the domain. We must 
impose one condition to select the Schwarzschild solution 
we wish to find, ie. to choose a mass. A simple way to 
fix this is that instead of solving the vv component of ([2| 
at 2; = 1, instead we replace it with a Dirichlet condition 

for T at Z = 1, so T|^^i = Tinner whcrC Tinner < tO 

ensure the domain pierces the horizon [45 . Consider the 
smooth metric. 



ds 



-(l-az)dt^ 



-dtdz 



VlO 



d^^ (12) 



z^ 

for constant a - note this is not Schwarzschild. As an 
example let us take this metric with a = 1.10 as the 
reference metric, and with a = 1.20 as the initial guess. 
The Newton method then converges to a solution with 
dzT\z=o = —1.06. Note that had one tried to find a 
black hole where the horizon was located at 2: > 1, for 
example taking an initial guess so that Tinner > then 
the method fails. The requirement of a smooth horizon 
is crucial to impose boundary conditions correctly. 

Two important points arise in this example. Firstly 
experimentally we find that for second order finite differ- 
ence the method fails, presumably as it doesn't impose 
smoothness of the functions to a sufficient degree. Cer- 
tainly for fourth order or above the method works very 
well, as it does for pseudo-spectral differencing. A second 
point is that if we had not imposed T = Tinner at z = 1 
but only the equations of motion there we would not ob- 
tain a locally unique solution. We emphasise that in this 
ingoing method the reference metric does not determine 
the moduli (in this case mass) of the solutions found. 



Appendix B: Details for inhomogeneous plasma 
flows and their duals 



For a given boundary metric deformation specified by 
the constants a and /3, the solution is determined by the 
parameters r, zq as discussed in the main text. We note 
that the global scaling = (t, p^y) ^ X of the bound- 
ary together with scaling of the parameters a ^ a and 
/3 A~^/3, and the boundary stress tensor components, 
{G4{Tab)) {G 4.{Tah)) relates solutions due to the 

conformal invariance of the boundary theory. In practice 
we choose zq = 1 and Zmax = 1.025, then vary r to obtain 
the required ingoing velocity Vq (or equivalently ingoing 
value of {T^P)/{T^^ + T^p)) which is scale invariant. In 
principle we would then use the above scaling to generate 
a solution with the required ingoing temperature Tq (or 
equivalently ingoing value of G4{Ttt)). However, for rea- 
sons that we do not understand (presumably related to 
the details of the way we fix the moduli of the solution) 
the value of Tq is actually the same to better than per- 
cent level for the various values of P we have explored, 
and so we have not needed to apply any scaling to the 
data presented here. 




FIG. 3. Coordinate positions of the horizon (solid lines) and 
ergosurface (dashed lines) for the flows in Fig. [l] The larger 
the value of /3 the deeper the horizon penetrates in the z- 
direction. Entropy is produced along the flow and hence the 
area density of the horizon is larger on the right end. 



In figure [sj we plot the position of the ergoregion (de- 
fined by T = 0), and the horizon solved from the o.d.e. 
gp^ d^hdyh = discussed in the main text. We note 
that both have a complicated dependence on x, although 
as expected this vanishes as X ^ ±1 where the metric 
becomes that of a homogeneous black brane, which in 
the coordinates defined by our reference metric will have 
ergoregion and horizon at constant z. Note also that the 
position of the horizon lies entirely within our domain 
< z < Zmax = 1.025 for all the solutions presented 
here. Since the metric functions are smooth at the posi- 
tion of the horizon, the future horizon is regular. 



We extract the vev of the CFT stress tensor from the 
bulk solution using the standard holographic renormali- 
sation prescription [41 . We proceed by first computing 
the boundary stress tensor in Feffer man- Graham coordi- 
nates. Then we construct, in a near boundary expansion, 
the change of coordinates from Feffer man- Graham coor- 
dinates to our working coordinates ([6| by requiring that 
= order by order. Defining, 



167r Ga 



Tu T, 



pt 



( T,, T, 







pp 







T 

^yy 



h{p) u^ip) 

uz{p) hip) 
ssip) 



(13) 

then ^3, 63, 1^3, 53 can be obtained from the bulk solution. 
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For example, for ts one finds an equation, 



^3(3;) 



1 1 
' 3 3! 

/4[-4a'(-6 



diT\ 



7/) + - 12xa(2))] (101s, + 5s3r) % 



(2) 



576cr 

-4xa'){387 5^ + 47 53,) 



1024 (j2 

f (cr')^(5981Sr + 993S3r) 



32768 

(14) 

where a is defined in ([5| and we give all functions 
in terms of the compact coordinate x (rather than p), 
f{x) = 1 — and s^r = sinh/cr. However, since the 
gauge condition = relates the derivatives of the var- 
ious metric functions, there are 3 other ways to extract 
ts from three z derivatives of the other metric functions 
than T. For a continuum solution these must all give the 
same answer, and we use this to test the accuracy of our 
stress tensor determination shortly. The other 63,113,53 
can similarly be extracted from three z derivatives of the 
various metric functions, and again there may be multi- 
ple ways to do this which are equivalent on Einstein solu- 
tions. The equations of motion also imply that 63, 113, S3 
are locally related to ts as the stress tensor is traceless 
and conserved, which again we check shortly. 

In figure [4] we display the vev of the T** component 
of the holographic stress tensor and compare it with the 
same component of the viscous hydrodynamic stress ten- 
sor, (16), for varying (3 in analogy with figure [l] in the 
main text. We see the same agreement with viscous hy- 
drodynamics at small /3, with strong deviations from it 
for P ~ 0(1). We note that (T^^) is well behaved even 
for the P = 2 fiow which has a superluminal region. 

In the figures above we have compared the holographic 
plasma behaviour extracted from the dual black holes to 
the viscous hydrodynamics predicted by the fluid/gravity 
correspondence. This gives a good approximation for 
\P/To\ <C 1. We now give details about these hydrody- 
namic fiuid fiows. Recall that the viscous fluid approx- 
imation to the plasma flow from fluid/gravity is deter- 
mined by the fluid stress tensor [46]. 



(15) 



-2 Q^t) P-P^^ (^V(,ix,) - ^^fi V^^/e) + O {\/'u) 

for a 1+2-dimensional boundary metric ds'^ = 
g^^^dx^dx^^ temperature T, 3- velocity with v? = —1 

and with Pah = UaU^ + g^^^ . The first term is that of an 
ideal fluid, and the latter is due to shear viscosity. For 
our flows, we take the 3-velocity in the p direction, so 
u'^ = (7, 7'u,0), with V the velocity and 7"^ = 1 — . 
For small but non- vanishing /3 the viscous term will gen- 
erate entropy and the fluid deviates from ideal behaviour. 




FIG. 4. Plot of (T") normalised by its ingoing value Tq vs. 
p for the holographic stress tensor (solid lines) and the stress 
tensor of viscous hydrodynamics (dashed lines). Top: Flows 
as in figure [1] with /3 = 0.2,0.3,0.5,0.7. Bottom: flows with 
P = 1., 1.5,2. For P 0(1) the stress tensor exhibits 0(1) 
features at small scales (compared to the length scale set by 
the temperature) and hydrodynamics no longer provides a 
valid description of the flow. 



The equation of motion for the fluid is given by conser- 
vation of this stress tensor. For our boundary metric 
there are two non-trivial components. One immediately 
yields G4{Ttx) cx l/^/o"- Combined with the second, one 
obtains flrst order o.d.e.s for the fluid velocity and tem- 
perature, as discussed in [30l [31]. It is these we have 
solved in order to compare to our numerical bulk solu- 
tions. In flgure[5]we display the behaviour of the temper- 
ature for the flows obtained from viscous hydrodynamics 
which we compared to the stress tensor from the gravity 
dual in flgures [l] and [4] ( we note that the velocity of these 
hydrodynamic flows are already shown in flgure [1]). For 
comparison with the viscous hydrodynamics we also show 
the ideal hydrodynamics solutions for the same ingoing 
fluid data. We see that as expected, for small f3 these 
closely agree, but for /3 ~ 0(1) the viscous behaviour 
departs strongly from the ideal behaviour, and likewise 
as we have seen in flgures [l] and [4j deviates from the full 
plasma behaviour as deduced from the gravity. 



Appendix C: Numerical errors and metric functions 

For the results presented in this paper we have dis- 
cretized the harmonic Einstein equation using sixth order 
flnite differencing, taking a uniform grid with lattice 
points in the z direction, and A^^^ = 4:Nz lattice points 
in X. We note that especially for the small /3 solutions 
where there are sharp gradients in the function a near 
the boundaries of the domain it is important to have suf- 
flcient resolution in the x direction to obtain accurate re- 
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FIG. 5. Temperature, T, plotted normalised by the in- 
going temperature To vs. p from viscous fluid/gravity hy- 
drodynamics (solid lines) and for ideal fluid/gravity hydro- 
dynamics (dashed lines) used to compare to the holographic 
plasma flows in figures [l] and [4] We emphasise that these 
only approximate the holographic plasma flow well for small 
13. Top: flows with /3 = 0.2,0.3,0.5,0.7. Bottom: flows with 
/3 = l.,1.5.,2. We see for small /3 agreement between the 
viscous and ideal hydrodynamics since there is little entropy 
generation. For larger /S ^ 0{1) viscosity becomes important 
and the behaviours strongly differ. 



FIG. 6. Convergence plots for the solution with f3 — 1. We 
show 8i ('o'), 82 ('x') and 83 ('+')? which measure the maxi- 
mum error in the Ricci scalar, Ricci tensor and the magnitude 
of respectively, as a function of the number of grid 

points in the ^-direction, Nz. The resolution in x is given as 
Nx — 4:Nz. We see linear convergence in this log- log plot. 
For 81 and 82 we find the slope is ^ 4 indicating fourth order 
convergence, whilst for 83 the slope is ^ 6. Other solutions 
exhibit the same convergence behaviour. 



suits. We have used resolutions up to A^^ x A^^; = 70 x 280 
points. Typically we begin by finding solutions at lower 
resolutions, and then use these as initial data to find the 
higher resolution solutions. 

We characterise the numerical error in our solutions by 
computing the error in solving the Einstein equations as, 



£1 = max 

0<z<H{x) 

£2 = max 

0<z<H{x) 



1 



12 



36 



(16) 



£s= max Jm^\ 

0<z<H{x) V 



where each is a scalar quantity which should vanish in 
the continuum for a solution, and is maximised for the 
solution in the exterior of the horizon. We maximise only 
over the region exterior to the horizon to obtain a well 
defined geometric quantity. We note that similar results 
are obtained when maximising over the entire domain. 
Typical results of convergence tests are shown for an in- 
termediate value of /3 = 1 in figure [6j For the maximum 
resolutions used, we see that the maximum local error in 
the solution is better than 10~^, as stated in the text. In 
addition, f 3 ^ in the continuum limit, which indicates 
that our solutions are not Ricci solitons. Similar results 
are obtained for the other values of P (including j3 = 2) 
discussed in this paper. 

We note that whilst we have used sixth order finite 



differencing, the slope of these curves against logTV^ is 
between ~ 4 — 6 depending on the quantity. We would 
naively expect ~ 6 for smooth solutions. We believe the 
observed lack of smoothness is not physical but due to 
our coordinate choice near the boundaries x = ±1. Since 
these boundaries are regular singular points of the p.d.e.s 
it maybe that there are {x ^ l)^log|x ^ 1| behaviours 
in the expansions of the metric functions for our gauge 
choice, where p is some power presumably with p > 4. 
We emphasise that we require only second derivatives 
to be defined for a solution to the Einstein equations, 
and the convergence we see is certainly much better than 
that, indicating the metric functions are better than 
in smoothness. As we show later, explicit calculation 
of the various two derivatives of metric functions gives 
well behaved results, again confirming better than 
smoothness. However, the apparent lack of smooth- 
ness does lead to poor convergence results when using 
pseudo-spectral differencing, hence our use of finite dif- 
ference. An obvious future direction is to improve the 
coordinate choice. 

We monitor the errors in the extraction of the stress 
tensor. As discussed in appendix B the components of 
the stress tensor t^{x)^ usix), ss{x) and bs{x) may be 
extracted from the 7 metric functions in different ways 
which should agree in the continuum. In figure [7| we 
display the fractional error in the different ways of ex- 
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FIG. 7. The function ts may be extracted from the metric 
functions in 4 independent ways which are equivalent in the 
continuum. Here we plot the maximum value of the three 
fractional errors, A^*-* characterising the deviation between 
these as a function of the number of grid points in the z 
direction, A^^, in a log- log scale for the — \ solution. We 
observe first order (or slightly better) convergence, which is 
consistent with the overall observed 4^^ order convergence as 
the quantity ts involves 3 derivatives of the metric functions. 
We obtain the same results for other values of /3. 



tracting ^3: 



A*^*^ = max 

-\<x<\ 



,(0) 



1,2,3, (17) 



where tl^^ denotes the value of ts obtained from d'^T\ 



as given in equation (|14[) , and 4*^ correspond to the other 



values obtained from the remaining independent combi- 
nations of metric functions. As this figure shows, A^*^ 
is consistent with vanishing in the continuum limit with 
a slope ~ 1 in a log- log plot. This is the expected be- 
haviour; from the equations of motion we have seen that 
we have fourth order convergence and the calculation of 
ts involves taking three derivatives of the metric func- 
tions. Therefore we expect the error in this quantity 
should exhibit approximately first order convergence. We 
see for our highest resolution data that the maximum 
fractional error is less than percent level. We obtain anal- 
ogous results for other components of the stress tensor 
which may be extracted in multiple ways, which we note 
includes the test of tracelessness of the stress tensor. 

Next we consider the error in the two non-trivial com- 
ponents of the conservation equation of the stress tensor: 



C2 



max 

-1<X<1 



max 

-1<X<1 



Uo 



Us 



_3 
^3 



2 cr 
/ 



1 



263 



(18) 



where each quantity should vanish in the continuum 
limit. In figure [S] we display Ci and C2 as a function 
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FIG. 8. Convergence plots for Ci,2, the maximum errors in 
the two non-trivial components of the conservation equation 
for the stress tensor for the /3 = 1 solution. The apparent 
linear convergence in this log-log plot against Nz has slope 
^ 4 for Ci, and ^ 1 for C2. These are consistent (or better) 
than expected, given the overall 4^^ order convergence, and 
that the conservation requires three derivatives of the metric 
functions in z and one in x. 



of the number of grid points in the z direction, A^^, again 
in a log- log plot. As this figure shows, the error in Ci 
exhibits almost fourth order convergence, which is bet- 
ter that one might have naively expected. On the other 
hand, the convergence in C2 is slightly better than first 
order, which is consistent with behaviour of the error in 
extracting and 63 as discussed above. We obtain simi- 
lar convergence results for the other j3 studied, including 
p = 2. 

In summary, our analysis of the numerical errors show 
that the solutions we present are of high quality, the 
maximum fractional error in the solutions being better 
than ^ 10~^. The finite differencing method we have 
implemented gives convergence to the continuum limit 
consistent with fourth order scaling. Our extraction of 
the components of the stress tensor exhibits the expected 
convergence, and we may estimate that the maximum er- 
ror in these components is better than 1%. 

We now turn to the metric functions. For concrete- 
ness, in figure [9] we show T and S over the domain for 
the P = 1 solution, and note that these are coordinate 
scalars with respect to z and x coordinate transforma- 
tions. In figure 10 we show the functions d^T and d^S to 
illustrate that the metric is better than and also that 
derivatives of the metric functions vanish as expected at 
X ^ zbl in the coordinate system defined by our reference 
metric (since the bulk solution becomes a homogeneous 
black brane there). The other metric functions show the 
same behaviour as those we show here. Likewise, tak- 
ing other two derivative combinations of these we obtain 
analogously well behaved functions. We obtain similarly 
well behaved metric functions for all other values of (3 
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presented here, including (3 = 2. 



FIG. 9. Metric functions T{z, x) and S{z, x) for the solution 
with (3—1. These and the other remaining metric functions 
V, B, F,U, A, are well behaved everywhere in our domain, in- 
cluding the region inside the horizon. 




FIG. 10. d^T and d^S for the /3 = 1 solution. As expected 
these are largest where a(x) changes most rapidly. We em- 
phasise that these and other two derivatives of the various 
metric functions are well behaved over the domain and van- 
ish at the asymptotic ends x ^ ±1 where the flows and dual 
black brane become homogeneous. 



Finally we plot the Weyl curvature as characterised by 



the scalar C^ypaC^^^^ over our domain in figure 11 and 



note that it is smooth, with no indication of any singular 
behaviour over the domain, again indicating the metric 
functions are better than smooth. 
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FIG. 11. C^^pcxC^^^" for the /3 = 1 solution. This function 
vanishes at the boundary of AdS and weh behaved elsewhere, 
indicating the absence of singularities. At x ^ ±1 where the 
flow and dual black brane are homogeneous the Weyl tensor 
also becomes homogeneous, with its x-derivative vanishing. 



